%% Characteristic Conditions
function [fields1,fields2] = chrt_bdry(fields1,fields2,model1,model2,layer1,layer2)

% fields1: the upper layer
nx1 = model1.nx;
% ny1 = model1.ny;
% dx1 = model1.dx;
% dy1 = model1.dy;
% dt1 = model1.dt;

lambda1 = layer1.lambda;
mu1 = layer1.mu;
% rho1= layer1.rho;
Zp1 = layer1.Zp;
Zs1 = layer1.Zs;
% no_of_lyr=1 means top
% no_of_lyr=n means bottom
% no_of_lyr1 = layer1.number;

vx1 = fields1.vx;
vy1 = fields1.vy;
sxx1= fields1.sxx;
sxy1= fields1.sxy;
syy1= fields1.syy;

% fields2: the lower layer
nx2 = model2.nx;
% ny2 = model2.ny;
% dx2 = model2.dx;
% dy2 = model2.dy;
% dt2 = model2.dt;

lambda2 = layer2.lambda;
mu2 = layer2.mu;
% rho2= layer2.rho;
Zp2 = layer2.Zp;
Zs2 = layer2.Zs;
% no_of_lyr=1 means top
% no_of_lyr=n means bottom
% no_of_lyr2 = layer2.number;

vx2 = fields2.vx;
vy2 = fields2.vy;
sxx2= fields2.sxx;
sxy2= fields2.sxy;
syy2= fields2.syy;

% normalization of nx1_axis & nx2_axis
nml_x1 = linspace(0,1,nx1);
nml_x2 = linspace(0,1,nx2);

% Always: nx1 > nx2
interp_vx2 = interp1(nml_x2,vx2(:,3),nml_x1,'pchip')';
% disp(size(interp_vx2));
interp_vy2 = interp1(nml_x2,vy2(:,3),nml_x1,'pchip')';
% interp_sxx2= interp1(nml_x2,sxx2(:,2),nml_x1,'pchip')';
interp_sxy2= interp1(nml_x2,sxy2(:,3),nml_x1,'pchip')';
interp_syy2= interp1(nml_x2,syy2(:,3),nml_x1,'pchip')';

vx1(:,end-2) = (1/(Zs1 + Zs2))*(Zs1*vx1(:,end-2)+Zs2*interp_vx2-sxy1(:,end-2)+interp_sxy2);
vy1(:,end-2) = (1/(Zp1 + Zp2))*(Zp2*interp_vy2+Zp1*vy1(:,end-2)-syy1(:,end-2)+interp_syy2);
old_syy1 = syy1(:,end-2);
syy1(:,end-2) = (Zp1*Zp2/(Zp1+Zp2))*(interp_vy2-vy1(:,end-2)+syy1(:,end-2)/Zp1+interp_syy2/Zp2);
sxy1(:,end-2) = (Zs1*Zs2/(Zs1+Zs2))*(interp_vx2-vx1(:,end-2)+sxy1(:,end-2)/Zs1+interp_sxy2/Zs2);
sxx1(:,end-2) = sxx1(:,end-2) + ((lambda1)/(lambda1 + 2*mu1))*(syy1(:,end-2)-old_syy1);

% downsampling
interp_vx1 = interp1(nml_x1,vx1(:,end-2),nml_x2,'pchip')';
interp_vy1 = interp1(nml_x1,vy1(:,end-2),nml_x2,'pchip')';
% interp_sxx1= interp1(nml_x1,sxx1(:,end-1),nml_x2,'pchip')';
interp_sxy1= interp1(nml_x1,sxy1(:,end-2),nml_x2,'pchip')';
interp_syy1= interp1(nml_x1,syy1(:,end-2),nml_x2,'pchip')';

vx2(:,3) = interp_vx1;
vy2(:,3) = interp_vy1;
old_syy2 = syy2(:,3);
syy2(:,3) = interp_syy1;
sxy2(:,3) = interp_sxy1;
sxx2(:,3) = sxx2(:,3) + ((lambda2)/(lambda2 + 2*mu2))*(syy2(:,3)-old_syy2);

fields1.vx = vx1;
fields1.vy = vy1;
fields1.sxx= sxx1;
fields1.sxy= sxy1;
fields1.syy= syy1;

fields2.vx = vx2;
fields2.vy = vy2;
fields2.sxx= sxx2;
fields2.sxy= sxy2;
fields2.syy= syy2;
end